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Abstract 

We analyze the dynamics of atom-laser interactions for atoms having multiple, closely spaced, 
excited-state hyperfine manifolds. The system is treated fully quantum mechanically, including the 
atom's center-of-mass degree of freedom, and motion is described in a polarization gradient field 
created by a three-dimensional laser configuration. We develop the master equation describing this 
system, and then specialize it to the low-intensity limit by adiabatically eliminating the excited 
states. We show how this master equation can be simulated using the Monte Carlo wave function 
technique, and we provide details on implementation of this procedure. Monte Carlo calculations 
of steady state atomic momentum distributions for two fermionic alkaline earth isotopes, ^^Mg and 
^'^Sr, interacting with a three-dimensional lin-_L-lin laser configuration are presented, providing 
estimates of experimentally achievable laser-cooling temperatures. 

PACS numbers: 42.50.Vk, 32.80.-t 
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I. INTRODUCTION 



The complex behavior that occurs when a multilevel atom interacts with polarization- 
gradient fields has been of interest for some time now. Sub-Doppler cooling ^ occurs 
because of elaborate optical-pumping processes produced by laser light in atoms with sub- 
level structure, as seen, for example, in the lin-±-lin and the (T+-cr_ laser configurations. The 
semiclassical understanding of these interactions more dimensions 

has led to a reasonably good qualitative understanding of the underlying mechanisms. Semi- 
classical analysis has even in some cases provided quantitative predictions of sub-Doppler 
laser cooling temperatures measured in experiments 

However, the most direct route to a quantitative understanding of atom-laser interactions 
is via a fully quantized master equation for the atom, in which the center-of-mass (CM) 
motion of the atom is taken into account quantum mechanically. This allows behavior at 
low laser intensities and low atomic velocities, the regime laser cooling strives to reach, to be 
described correctly. The drawback of solving such a master equation, however, is the large 
number of basis states required for the calculation, due to the additional momentum states. 
This problem becomes especially pronounced when attempting to model three-dimensional 
(3D) systems, where the state space grows as the cube of the number of one-dimensional 
momentum states needed. 

The Monte Carlo wave-function (MCWF) technique, introduced in the early 1990's has 
allowed significant progress to be made on the subject of atom-photon interactions in 3D 
as well as lower-dimensional calculations. The MCWF technique is a simulation procedure 
for the master equation that involves propagation of single stochastic wave functions, rather 
than density operators, with random processes occurring at random intervals due to inter- 
actions with the photon field that cause spontaneous emission. It has been shown that this 
method is equivalent to the master equation in the limit of a large number of independent 
stochastic wave functions j^. The MCWF technique has been successfully utilized to cal- 
culate 3D sub-Doppler laser cooling temperatures for atoms with Zeeman degeneracy in the 
ground and excited states j^. 

The majority of the research done on laser cooling has involved essentially two-level 
systems, consisting of a ground state and an excited state, which may or may not contain 
degenerate sublevels. However, some investigations have explored atomic systems in which 
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multiple distinct excited states come into play. In particular, the use of bichromatic laser 

(seeRefs. 

for example). 

This paper focuses primarily on monochromatic laser cooling for atoms with multiple 
closely spaced hyperfine excited-state manifolds. Figure ^ provides a graphical illustration of 
this type of atomic configuration. This situation is of importance, for example, in alkaline- 
earth atoms with nonzero nuclear magnetic moment. If the excited state manifolds are 
spaced in energy on the order or smaller than the excited state linewidth 7, coherences 
between these manifolds become nonnegligible, and can have a significant effect on the optical 
pumping processes required for sub-Doppler cooling and on the dynamics of the atom-photon 
interaction. Sub-Doppler laser cooling was experimentally identified in fermionic '^^Sr jisj ], 
despite significant spectral overlap in the excited state. At the time, it was hypothesized 
that the large ground-state degeneracy in ^"^Sr (due to the large nuclear spin / = 9/2) 
was somehow able to overcome the decrease in cooling due to the spectral overlap. Other 
systems with s pec tral overlap in the excited state are Q, ^Li Q, and the fermionic 
isotopes of Yb In ®^Rb, the effects of excited-state spectral overlap on the effectiveness 
of velocity-selective coherent population trapping have been explored, both experimentally 
and theoretically Our goal in the paper is to provide a detailed discussion of the 

theoretical techniques required to model such systems realistically. In a future publication, 
we plan to present comprehensive laser-cooling predictions for a variety of atoms. 

The structure of this paper is as follows. In Section II, we develop the master equation for 
a laser-driven atom with multiple excited-state manifolds, and then specialize this equation 
to the low- intensity limit. In Section III, we introduce the MCWF technique and apply it 
to this low-intensity master equation. In Section V, we perform full Monte Carlo master- 
equation simulations for ^^Mg and ^^Sr atoms in a 3D lin-±-lin laser configuration as an 
example of using this technique determine expected temperatures for these atoms in a laser 
cooling experiment. In Section VI, we conclude. 

II. MASTER EQUATION IN THE LOW-INTENSITY LIMIT 

In this section we develop the master equation describing a multilevel atom interacting 
with a coherent laser field and coupled to a vacuum photon field. It is this equation, with 
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FIG. 1: Energy level diagram of an atom with multiple hyperfine manifolds. If the energy spacing of 
the excited-state manifolds are of the order or smaller than the natural linewidth of the transition, 
the usual sub-Doppler cooling transition (Fg <^ F^. = Fg + 1) is not isolated and the other manifolds 
must be taken into account. 

quantized atomic CM, that will provide an accurate description of atom-photon dynamics, 
and this master equation will provide the basis for the Monte Carlo simulations that will be 
discussed later. 

The full Hamiltonian for the atom-laser system plus the radiation field is 

H = Ha + Hr + Va^l + Va-r, (1) 

where Ha = ^i^i + |^ bare atomic Hamiltonian, is the vacuum radiation 

field Hamiltonian, and Va-l and Va-r are the atom-laser and atom-radiation field coupling 
terms, respectively. In the atomic Hamiltonian, Pi is a projection operator onto the i-th 
internal excited-state manifold, fiWi is the energy of the i-th excited-state manifold relative 
to the ground-state manifold, P is the atomic CM momentum operator, m is the atomic 
mass, and the sum runs over all excited-state manifolds. We have assumed in Eq. (HI) that 

n 

the effects of atom-laser and atom-radiation- field coupling are independent |20||. 

We can view Eq. in terms of system-reservoir interactions. The system consists of the 
atom, the laser, and their interaction. The system Hamiltonian is 

Hs = Ha + Va-l. (2) 

The reservoir is the vacuum radiation field, having many more modes than the system. With 
the Markov approximation, along with a few other approximations, the master equation is 



then given by 

&^^[a,Hs]+jC,p[a]. (3) 

The operator a is the system reduced density operator element, i.e., the reservoir degrees 
of freedom have been traced over, a = Tvup. The remaining term, £sp[(T], encompasses the 
interaction between the atom and the vacuum photon field, and provides for the phenomenon 
of spontaneous emission. 

The relaxation operator due to spontaneous emission, which we derive in detail in the 
Appendix, is given by 

-pL f£nJ2Y.[(^- A^^) ^^■''e-'''-^(e* • A(^))(7 + a{e ■ A« V'^'''e-^'^-^(e* • A^^))] , 

^ e±k i,j 

(4) 

where A^*)^ and A*^') are vector raising and lowering operators, respectively, between the 
ground state and the ith excited state, R is the atomic CM position, k is the direction of 
the photon emitted in the relaxation process, and 7 is the decay rate of the exited states. 
The integral is performed over solid angle in the vector k and the sum over e ± k refers to the 
two polarization directions perpendicular to k. Note that here and throughout this paper, 
we assume that each of the excited-state hyperfine manifolds has the same lifetime r = 7"^. 
Expanding these vector operators in a basis of spherical unit vectors, e±i = ^{x ± iy)/y/2 
and eo = z, we have 

q=0,±l 

The spherical components of the vector operators are 

<^F„F.„M.,Me„w\J9lF9Mg){JeIFeM (6) 
E <^F,,F.,M.,Me^,J,,J.,l\JeIFeM,,){J,IF,M,\, (7) 

Mg,Me. 

where 

«F„F„M„M,,.„..,/ = (-l)We,+M,+..+7^(2F, + 1)(2F, + 1)(2J, + 1) 
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Eq. (0} is written in a way that makes explicit that it is in Lindblad form [2l|, |22 
As we will see later, it is important for the relaxation operator to be of this form in order 
to make use of the MCWF technique. Because the complex exponentials in the second line 
cancel each other, the remaining integral over solid angle can be evaluated, whereby Eq. (@)) 
can be equivalently written as |8| 

e±k i,j i 

We will now examine that atom-laser interaction term, which is given in the electric-dipole 
approximation by 

VA-L(R,t) = -D-Ei(R,t), (10) 

where Ei(R, t) is the electric field of the laser and D is the electric dipole operator. As usual, 
we treat the laser as a classical field, since it is a densely populated mode of the electric 
field. We can write the laser electric field in terms of its positive and negative frequency 
components, El(R, t) = E^^''(R)e^*'^* + c.c, and then expand into spherical components, 

Ei:^(R) = Y E (-i)X(R)^-., (11) 

g=0,±l 

where Eq is the electric-field amplitude and aq(R) are the expansion coefficients. Making 
the rotating-wave approximation, so that 

Va-l{R, t) = -DW ■ Ei+^(R)e-^"* - D(-) ■ Ei^"^(R)e^'^*, (12) 

where D(+) = ^. Pe^'DPg and D^") = PgDPe,, we find 

Va^l = -| E + H.c. (13) 

i 

In the previous equation we have defined the atom-laser raising operator, 

viiR) = j2 (14) 

g=0,±l 

and lowering operator, 

P.(R) = al{R)A^\ (15) 

q=0,±l 

and introduced the "invariant" Rabi frequency, 

Q = Eq^ J^W^W ^ (16) 
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where {Je\\D\\Jg) is the reduced dipole matrix element between the ground and excited 
states. This form of a Rabi frequency, defined in terms of the reduced matrix element 
between the J = Jg ground state and the J = excited state, is convenient because, in 
general, Rabi frequencies for transitions to different excited-state manifolds will not be the 
same. 

Next, we observe that the second term in Eq. Q is comprised of excited-state projection 
operators both pre- and post-multiplying the system density operator. Thus, it is clear that 
this term can be absorbed into the free-evolution commutator term in Eq. Q, allowing the 
master equation to be equivalently described by Hamiltonian evolution determined by an 
effective Hamiltonian ifgfr, plus a term which is commonly called a jump term, and which 
cannot be written in the form of a commutator with the system density operator. We thus 
have, 

a = [H,,a - aHl,) + g / e-'^-^(e* ■ A«)a(e • A^^)^^''-^, (17) 

where the effective Hamiltonian ifcfr is given by 

i 

where Va-l is as given in Eq. In obtaining Eqs. (fT7j) and p8|l . we have made the usual 
rotating-frame transformation, which removes the free-evolution atomic Bohr frequencies 
from the problem. The more relevant frequencies are instead the laser detunings 6i = uj — Ui 
from the ith excited-state hyperfine manifold. The master equation given in Eq. ()17|) is 
fully general, but has been written in a form that will facilitate setting up a stochastic wave 
function simulation using the MCWF technique described later. 

We would like to now specialize the master equation just discussed to the limit of low 
laser intensity. Specifically, this limit is valid when the saturation parameter for the atom 
in the ith excited-state hyperfine manifold, 

- WTWW' ^^^^ 

is small, which occurs when the laser intensity is small or the laser detuning from the atomic 
transition is large. In this limit, the excited states are said to adiabatically follow the ground 
states. The excited states can then be eliminated from the equations of motion, resulting in 
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a master equation in terms of only the ground-state sub-density-matrix, 



= Pg^Ps- (20) 
In this hmit, the master equation becomes (see section 8.3.3 of Ref. ^24] ) 

a,, = ~ (W,, - a,,hl^) + fd'QYl ■ k))a,,(e ■ B« V, k)). (21) 

The new effective Hamihonian is given by 

^eff = ^ + E (^^ - (22) 

i 

The new decay raising and lowering operators are given by 

Bf{K, k) = y^A»^e^'^-^P«(R), (23) 

and 

5«(R, k) = ./5240e-k.RpWt^j^^_ ^24) 
^ V 87r 

Note that this new lowering (raising) operator contains two components: a raising (lowering) 
operator D(*)^(R) ('D*^*)(R)) between the ground state and the ith excited-state manifold due 
to the atom-laser interaction, and a lowering (raising) operator A^q^"^ e^^'^ {Aq^ e'"^^'^) of type 
q corresponding to coupling with the reservoir photon field via a photon with polarization 
q. Thus, the jump operator in the low-intensity equations describes a transition cycle of 
the atom involving coupling to both the laser and the reservoir photon field. Note also that 
this new operator and the effective-Hamiltonian term in the equation of motion are both 
proportional to the saturation parameter Sj, the perturbation parameter. 



III. THE MONTE CARLO WAVE-FUNCTION TECHNIQUE 

ISa 1271 |28|, 



The MCWF M 



29| technique is a means of interpreting a system- 



reservoir master equation — which describes the evolution of a density operator for a system 
interacting with a large external reservoir — as the evolution of an ensemble of individual 
wave functions, each undergoing random quantum jumps. The free evolution of the stochastic 
wave functions is determined by the effective Hamiltonian that we found in the previous 
section. The nature of the quantum jumps is determined by the leftover term in the master 
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equation, which cannot be absorbed into the free-evolution commutator. The components 
of this leftover term are often called quantum-jump operators. 

In the following, we will deal primarily with the master equation in the low-intensity limit, 
as developed in the previous section, although the methods could just as easily be applied 
to the arbitrary-intensity master equation. The low- intensity limit, however, provides a 
reduction in the number of internal atomic states required in the calculation, and this will 
be beneficial for performing calculations later. Furthermore, since the lowest temperatures 
are achieved for low laser intensities, such a specialization does not hinder our ability to 
calculate lower bounds of temperature. 

Having already expressed the master equation in a form involving an effective Hamiltonian 
and a jump term in the previous section, the application of the MCWF technique is rather 
straightforward along the lines developed in the literature (see, in particular, Ref. For a 
single stochastic wave function, the procedure is as follows. First, set the wave function to an 
initial value. Then, numerically propagate the wave function for a time step 5t according to 
the effective Hamiltonian ifefr only, from an initial value \ip{t)) to a final value \ip'^^\t + 5t)), 

|^«(t + 5t)) = (l-^)|^(t)). (25) 

Restrictions on the size of 6t are given such that the first-order truncation of the time- 
evolution operator in Eq. ()25|1 is approximately valid. We note that Hes is non-Hermitian 
by construction, as a result of absorbing parts of the relaxation operator into the original 
(Hermitian) bare system Hamiltonian. Because of this, propagation with Hcs will not con- 
serve the norm of the wave function when propagated to {ip^^^t + St)). The time step 6t of 
the propagation must be chosen so that 5p ^ 1 in the inner product, 

(V^(^) {t + 6t) \^^^\t + 6t)) = 1 - 6p. (26) 

The quantity 6p is the loss of norm resulting from propagating with i^gfr for a time step 6t, 
and is found to be 

5p = 5t (V^(t)l 5^BW\R,k) ■ B«(R,k)|V^(t)) 

i 

= 6t{m\Y. E Bf{K,k)Bj;\K,k)\m) (27) 

i q=0,±l 

i g=0,±l 
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The total loss of norm has been decomposed into individual elements each corresponding to 
a particular type of interaction with the reservoir (i.e., the g- value of the interaction, or the 
excited state i involved). These individual contributions are given by 



We see that the loss of norm due to a given type of interaction with the reservoir is determined 
by the quantum-mechanical expectation value of the product of jump operators of this type 
of interaction. The loss of norm 6p can also be interpreted as the probability for a quantum 
jump to occur. 

After the wave function has been propagated as described above, and the values of Spi^q 
calculated, it must then be determined whether or not a quantum jump occurred. This is 
achieved by generating a pseudo-random number on a computer and comparing it to the 
value of the total jump probability 6p. If the random number is less than 6p, a quantum 
jump occurred, and if it is greater, no quantum jump occurred. If a quantum jump does 
occur, the type of quantum jump must also be calculated by comparing the random number 
with the individual sub-probabilities 6pi^q in the same manner. 

If a quantum jump of type q, i occurs, we must apply the quantum jump lowering operator 
Bq^\'R, k) to the wave function from the beginning of the time step. 



The square-root factor in front of the lowering operator is necessary for renormalization. If 
no quantum jump occurs, then we simply renormalize the wave function. 

The resulting wave function is then used as the starting point for propagation over the 
next time step, and the procedure is repeated. 

A good approximation of the true system density matrix is achieved by combining the 
trajectories of a number of independently propagated stochastic wave functions, each tra- 
jectory having a unique sequence of pseudo-random numbers. (A thorough discussion of 
the statistical issues involved with the MCWF technique can be found in Ref. Q.) Once 
a suitable ensemble of stochastic wave function trajectories has been obtained, an estimate 
of the true expectation value of an operator is found by taking the ensemble average of 
the expectation value of that operator with respect to the stochastic wave functions. For 



6p,,q = 6t{m\B^^ 



(R,k)i?»(R,k)|^(t)). 



(28) 




(29) 
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example, an estimate of the average kinetic energy at a time t for a system for which 
independent stochastic wave functions have been calculated is given by 



where ipi{t) is the ith stochastic wave function, given at time t. 

Figure |21 demonstrates a simple example of the application of the MCWF technique, 
wherein the average kinetic energy is calculated for a two-level atom interacting with a 
one-dimensional standing-wave field. For this calculation, we have used a Rabi frequency 
of Vt = 7/2 and a detuning of 6 = —7/2, where 7 is the decay rate of the upper to the 
lower atomic state, and we have set 7 = 400-E'r., where Er = h^k"^ /2m is the recoil energy. 
The atomic kinetic energy, averaged over 500 stochastic wave functions each initialized to 
zero momentum, is plotted as a function of time, with error bars indicating the error in the 
ensemble average for a given time. The separation of the transient relaxation period from 
the steady-state is clear, the steady state regime being characterized by fluctuations in the 
average energy about a mean. This noise is due to the flnite number of stochastic wave 
functions being used, and if a greater number of wave functions were used, the amplitude of 
the fluctuations would be decreased. In the limit of an inflnite number of wave functions, 
the true density-matrix solution of the master equation would be obtained. An estimate 
of the steady-state kinetic energy is found by time-averaging the calculated data over the 
entire steady-state regime. Since this is a larger ensemble than the set of wave functions for 
a single time, the error of such an average will be smaller than the error bars shown in the 
flgure. 

IV. CALCULATIONS FOR ^^mG AND ^^SR 

The purpose of this section is to illustrate the application of the theory developed up 
to this point to a complicated system. We wish to quantitatively study the dynamics of 
particular atoms interacting with 3D polarization-gradient laser flelds. The balance of the 
frictional cooling forces along with the diffusion experienced by the atom due to spontaneous 
emission and its interaction with the laser leads to a steady-state momentum distribution 
that determines the temperature of a gas of such atoms. In particular, we will study here the 
cooling of the fermionic isotopes of two alkaline-earth atoms, ^^Mg (nuclear spin / = 5/2, 
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FIG. 2: (Color online) An example of a characteristic MCWF stochastic trajectory. Shown is the 
average of the atomic CM energy over 500 independent stochastic wave functions, as a function of 
time, for a two-level atom in a ID standing-wave laser field. The energy is given in units of the 
recoil energy Er = H'^k'^ /2m, and time is given in units of the inverse recoil frequency = h/E^. 
All wave functions are initialized in the ground state of the atom and localized in momentum space 
with zero momentum. The steady state, wherein the system fluctuates around an average value, 
is seen to be achieved after a transient relaxation period. Error bars indicate the variance in the 
data at each given time for the ensemble of 500 stochastic wave functions. An estimate of the 
steady-state atomic CM energy is obtained by performing a time-average over all wave functions 
for all times after the relaxation regime. The error bar of such an average will be smaller than the 
error bars in the figure, which apply only to the data for a given time. 

^So-^Pi width 7/27r = 81 MHz, hyperfine splittings Au;i3/27r = 46 MHz and Awaa/Svr = 
27 MHz, where we have assumed a hyperfine quadrupole parameter i? = [s^) and ^''Sr 
(/ = 9/2, ^Sq-^Pi width 7/27r = 32 MHz, hyperfine splittings Acji3/27r = 43 MHz and 
Aco'23/27r = -17 MHz). These atoms, having nonzero nuclear magnetic moment, have de- 
generate (assuming zero magnetic field) Zeeman sublevels. These sublevels allow for the 
mechanism of sub-Doppler cooling in an appropriate laser configuration. Both ^^Mg and 
^'^Sr exhibit significant excited-state spectral overlap, with Ac<Ji3/7 = 0.57, /S.uj2zh = 0.33, 
and Aci;i3/7 = 1-3, Aci;23/7 = -0.53, respectively. We consider the 3D lin-±-lin laser configu- 
ration, consisting of a pair of opposing beams along each cartesian axis, in which each beam 
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is linearly polarized orthogonal to its opposing beam. Furthermore, for this calculation, we 
set to zero the relative phases of the three sets of laser pairs. 

Having a nuclear spin of / = 5/2, the ^So state of ^^Mg results in a hyper fine ground 
state with 6 sublevels. Use of the low-intensity master equation given in Eq. (jl7p allows 
us to consider only these 6 internal states of the atom, since the excited states have been 

Q 

adiabatically eliminated in this regime. However, as noted in Ref. [9|, a momentum grid 
extending to 20hk in each direction with a spacing of hk would yield a density matrix 
with (6 X 41'^)^ ^ 2 X 10^^ elements. A direct solution of this master equation is not 
numerically feasible, even without considering the further increases in matrix size necessary 
to describe the master equation relaxation operator in Liouville space [s])]. On the other 
hand, the MCWF method only requires numerical propagation of individual wave functions, 
which would be represented by vectors with 6 x 41^ ^ 4 x 10^ elements. If the number of 
independent stochastic wave functions required to achieve satisfactory convergence for the 
calculation of a particular property of the system is not unreasonably large, the MCWF 
method provides a distinct advantage over a direct master-equation solution. 

We follow the procedure outlined in Section HI, working in the low-intensity limit in 
order to reduce the number of internal atomic states in the calculation, which increases the 
efficiency of calculation. Since laser cooling is most effective at low laser intensities, this 
turns out to be a useful regime in which to work, with the additional benefit that lower 
temperatures require a smaller number of atomic CM momentum states in the calculation. 
We must determine the effective Hamiltonian as given in Eq. ()22|1 and the jump operators as 
given in Eqs. and for each atom, and for the particular laser field being considered. 

We consider here the lin-±-lin laser configuration in 3D, with the relative phases of the 
beams set to zero. The positive-frequency component of the electric field is 




(31) 



9=0,±1 
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with spherical coefficients 

a+i(R) = (e-''^^ + e*^^ + ze''^^ + ze-''^^) , (32) 

a_i(R) = (e-^^^ + e*^^ - ze*^^ - ^e-''^^) , (33) 

ao(R) = e-^'^^ + e'^^. (34) 

With these coefficients, along with parameters appropriate to the particular atom under 
consideration, the atom-laser raising and lowering operators given in Eqs. (fT^ and ()15p can 
be constructed. With knowledge of the effective Hamiltonian and the raising and lowering 
operators, we can then proceed with the MCWF procedure as outlined. 

Our example entails propagating 20 stochastic wave functions each for three different 
values of the light-shift parameter, h\5z\s^/ {2E,^c) =10, 20, and 30, for both ^^Mg (/ = 5/2) 
and '^'^Sr (/ = 9/2). We consider only ^3 = —87. As in Figure |21 we calculate the stochastic 
trajectories of the ensemble average (i.e., averaged over the 20 wave functions) kinetic energy 
for each atom as a function of time. We continue this propagation until the transient regime 
has been passed for some time, and use the time average over the steady-state ensemble- 
average kinetic energy to provide an estimate of the total average kinetic energy and the final 
error. The results are shown in Fig. IHl along with the energies for atoms with an isolated 
cooling transition for comparison, Jg = Jg+^ with Jg =1, 2, 3, and 4, with detuning 5 = —87, 
as ffist calculated by Castin and M0lmer in Ref. 9|. From this cursory analysis, we can see 
that ^^Mg should exhibit a sharp rise in temperature with increasing laser intensity, while 
^^Sr will cool to sub-Doppler temperatures even for higher intensities, as has been noted 
experimentally jl^. 

Detailed calculations of this sort, for realistic atoms, are quite computationally expensive. 
For example, a single data point for the Mg and Sr calculations presented here required on 
the order of 200 hours wall time for a 20 processor parallel code, running on a cluster of 
2.4 GHz Intel Zeon processors. There remains work to be done improving the numerical 
efficiency of our initial codes. Our goal in this paper has been to present our method 
and some illustrative results; in a future publication we plan to expand upon these initial 
results using improved, faster codes and present comprehensive predictions of laser cooling 
temperatures for a variety of atoms. 
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FIG. 3: (Color online) Results for calculated ensemble-average energies (rms momentum squared) 
for ^^Mg and ^''Sr, as a function of the light-shift parameter h\63\s3/{2Ej-cc)- For comparison, also 
shown are the calculated energies for atoms with isolated transitions, Jg = Jg + 1, with Jg =1, 2, 
3, and 4, with detuning 6 = —57. See text for discussion. 

V. CONCLUSIONS 



In conclusion, we have provided a detailed description of the fully quantum-mechanical 
master equation that describes an atom with multiple internal internal structure interacting 
with a 3D polarization-gradient laser field. We have shown how the spontaneous-emission 
relaxation operator is generalized for atoms of this type. The MCWF technique has been 
applied to these equations of motion, providing a more efficient means of performing calcu- 
lations for these systems compared to a full solution of the master equation. A few example 
calculations have been presented to illustrate the application of this theory to atomic sys- 
tems interacting with laser configurations commonly used in experiments. After making 
improvements in the efficiency of our codes, we intend to expand upon this work in a future 
publication and provide a comprehensive survey of laser cooling calculations for atoms with 
multilevel internal structure. 
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APPENDIX: RELAXATION OPERATOR FOR AN ATOM WITH MULTIPLE 
EXCITED-STATE HYPERFINE MANIFOLDS 

In this appendix, we outline the major steps in deriving the spontaneous-emission re- 
laxation operator for an atom with multiple hyperfine excited-state manifolds. Detailed 
derivations of this sort, but including only a single excited state manifold, exist elsewhere 
in the literature (see, for example, Ref. [3^). Our intent here is to highlight the steps im- 
portant in generalizing the previous work to include coherences between other exited states. 
We will work within the framework of the theory of system-reservoir interactions and follow 
the notation of Ref. 2^ . 

The total Hamiltonian for an atom coupled to a vacuum radiation field is then given by 
H = Ha + + Va-r, where Ha and H^ are the atom and reservoir bare Hamiltonian, 
respectively, and Va-r is the atom-reservoir coupling and is given in the the electric-dipole 
approximation as 

Va-r = -D ■ E(R) = - J2 i-iyD^E.^iR). (A.l) 

(j=0,±l 

Here, D is the electric dipole operator for the atom, and E is the electric-field operator for 
the photon field, and we have expanded the interaction into its spherical components. To 
simplify the formalism, we will begin by ignoring the atomic CM and setting the position 
coordinate to be the origin, R = 0. At the end we will then generalize the equations to 
include the CM degree of freedom. 

In general, the total density operator p evolves according to the Liouville equation, p(t) 
i [p (t), H]. Making the usual assumptions involved in deriving the master equation j2 
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l32l |. we arrive at an equation of motion for the reduced density operator of the system 
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'\ poo 

a{t) = -[a{t),H^]-- rfr^(-l)'' 

In the previous equation, gq{T) is the two-time correlation function of the reservoir and 
is defined as gqi^r) = TTji[aREg{T)Eg{0)], where the variables with tildes are operators in 
the interaction representation, Eg{t) = e^^R^/f^EqC'^^'^^/^. We assume that the reservoir is 
initially a vacuum, so that Or = |0) (0|. From this we can see Qq^r) = K^|-^g|0)l^6~*'^''*) 
where the kets and bras refer to reservoir states. Note that 5'(t)* = g{—T). The correlation 
time of the reservoir tc is defined such that g{T) ^0 for r ^ re- 
in addition to the above approximations, we will also make the secular approximation, 
which requires that the equation of motion for each density-matrix element dij have only 
terms involving density-matrix elements a^i on the right-hand side such that \ujij —uJki\ <^ 7, 
where and where 7 is the order of magnitude of the system-reservoir coupling. 

In the following, we will consider a system with a ground state coupled to multiple excited 
states that are separated in energy of the order or smaller than 7. Thus, the ground-excited 
energy splitting \ujge\ ^ 7 will be a non-secular frequency, while oj^iej ~ 7 will be a secular 
frequency. 

The particular atomic system that we are considering consists of an ground state with 
electronic angular momentum J = Jg = and an excited state with J = Jg = 1. The 
electronic angular momentum is coupled to the nuclear spin quantum number J, resulting 
in a ground state with total angular momentum Fg = I, and three excited states with 
{Fe-} = {1 — 1,1,1+ 1}. These assumptions are made for concreteness, but we note that 
this derivation can be easily extended to arbitrary angular momentum schemes. It is useful 
to decompose the system density operator as illustrated in Fig|l[ 

where (Tij{t) = Pi(j{t)Pj; Pi is a projection operator onto the i-th hyperfine manifold. Pi = 
YIm' \ J^P^i) {JIF^i\'i and Mi is the substate label for the i-ih manifold. Two relations 
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FIG. 4: The partitioning of the density operator for an atom with multiple coupled excited-state 
manifolds, each potentially having multiple substates. 



that will be useful in the following are 

J2 (M,|Z},|M,J|M,)(MeJ=4 



^) {J9\\D\\Je) 
V2 + 1 ' 

^^(i) HJg\\D\\Je) 

" V2 Je + 1 



(A.4) 
(A.5) 
), and 



where Aq"^^ and Aq"^ are the atomic raising and lowering operators defined in Eq. 
where we have made use of symmetry properties of the three-J and six-J symbols 

We focus on the equation for the ground-state sub-density-operator <Jgg{t) in Eq. ()A.3|) . 
Beginning by taking matrix elements of Eq. ()A.2|1 between ground-state sublevel kets, we 
proceed as usual by eliminating terms that violate energy conservation (e.g., photon emis- 
sion coupled to atomic excitation), and absorbing interaction- induced energy shifts into the 
energies of the internal atomic levels. The resulting equation of motion is 

3 1/ T llr^M T\|2 



3eo(27r)3c3n 



2 Je + 1 ^ ^ 1 ' 1 

1 «.i 



7 



(A.6) 



We have assumed that the energy splittings between ground state and the various excited 
states are all approximately equal, and accordingly have defined ujq = uj^.. — Ug for i =1,2,3. 
Equivalently, we have assumed that the decay rate for all of the excited-state manifolds is 
approximately equal, and defined 7 = 7j^.^j for i =,1,2,3. Note that the double sum over 
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excited-state manifolds in Eq. ()A.6|) will clearly result in inter-manifold coherence effects in 
the equations of motion. 

Regarding the energy shifts of the internal atomic states that arise due to interaction 
with the reservoir states, it is important to mention a subtle feature not found in the sim- 
pler case of degenerate isolated excited states. Such energy shifts occur in the form of 



divergent principal-part integrals of virtual transition amplitudes l32| |. For degenerate 
isolated manifolds, these diverging terms can be shown to cancel each other in the equa- 
tions of motion. However, for the case of multiple, nondegenerate manifolds, these terms 
no longer cancel exactly, and pathological divergences related to reservoir-dressed internal 
atomic energy splittings remain. A thorough exploration of these terms is outside the scope 
of this paper, and for the present purposes, we ignore such diverging terms and absorb 
interaction-induced energy splittings into the defined energy levels of the atoms. 

Working in the same manner as for the ground-ground sub-density-operator, we can find 
the equations of motion for the excited-state sub-density-operators, 

= -za;,,^P,a(t)Pe, - | E E A^<l^^{t) + cr{t)4'^U?) P^,, (A.7) 

g k,l 

where we have added a trivial summation index that will be useful later when combining the 
various sub-density-matrix decay terms. Similarly, the equations of motion for the optical- 
coherence sub-density-operators are 

de,,(t) = -tU^^.P^^amj - i E E ^e.AW^4V(t)P„ (A.8) 

q k,l 

and &geM = ^lg{t). 

Using Eq. ()A.3|) . we can construct the equation of motion due to spontaneous emission 
for the full density operator. 



(A.9) 



Defining the spontaneous emission relaxation operator. 



we can write the equation of motion as 



^{t)='-[a{t),HA] + C,p[cT]. (A.ll) 
19 



Including the atomic CM dependence that we have been ignoring since the beginning 
amounts to adding an integral over momentum states in 3D that should have been included 
when we inserted atomic projection operators. With this addition, the full relaxation oper- 
ator takes the form shown in Eq. 1^. 
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